Image lens distortion correcting method

ABSTRACT

An image lens distortion correcting method having step (A) for printing an arbitrary image I 1  in a computer; step (B) for photographing the printed image I 1  with a camera having a lens distortion and obtaining a photographed image I 2  in the computer; step (C) for obtaining a parameter θ such that the image I 1  is artificially distorted with the parameter θ and an obtained distorted image I 1   ud  equals with the photographed image  2 ; and step (D) for using the obtained parameter θ to correct the image photographed by the camera. Thereby, (1) no correspondence point needs to be given, (2) no special tool or object is used, and (3) all points on the image are used, so that an internal parameter can automatically be obtained with high precision, and image lens distortion can be corrected with high precision.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to a correcting method of lens distortionof an image.

2. Description of Related Art

An image taken by a video camera or a digital camera causes distortion(distortion aberration) depending upon a camera lens. The distortionaberration is generally extreme with a wide angle, increases toward aperiphery of the image, and increases to 3% or more particularly with azoom lens. Moreover, since the image is deformed in a so-called barrelor bobbin shape, this distortion aberration disadvantageously adverselyaffects estimation in computer vision, measurement using the image, andthe like.

Parameters such as a distortion coefficient for determining a cameralens distortion are called a camera internal parameter, and cameralcalibration means that the camera internal parameter is obtained. Whenthe camera internal parameter can be obtained through cameracalibration, the image lens distortion can be corrected by an imageprocessing using the internal parameter.

A conventional camera calibrating method is disclosed in many researchpapers.

In the aforementioned conventional method, a relation between a certainpoint in a three-dimensional space and a point signal to the certainpoint (correspondence point) in the image is derived, and the internalparameter is estimated. There are two methods of obtaining thecorrespondence point as follows:

(1) A person designates the correspondence point with a mouse or thelike.

(2) A grid intersection point, a polygon or polyhedron apex, or anothercharacteristic point is automatically detected.

Since either one of these methods includes a measurement error, a largenumber of correspondence points have to be given in order to inhibit anerror influence. That is, the more the number of correspondence pointsis, the more precise the obtained internal parameter becomes. Therefore,correspondence of an enormous number of (several tens to severalhundreds) points has to be designated.

Moreover, these methods have the following problems.

In the method (1) in which the person designates the correspondencepoint with the mouse or the like, the error is large, operation ismonotonous, and much labor is necessary. Furthermore, fatigue increasesthe error during designation.

In the method (2) (characteristic point detecting method) ofautomatically detecting the characteristic point by the imageprocessing, the number of points detectable by the processing islimited. For example, when the grid intersection point is detected, thenumber of intersection points is only 88 in total even with 8longitudinal lines and 11 lateral lines. It is said that at least 200points or more are necessary for obtaining the parameter with goodprecision, and it is necessary to prepare the corresponding grid,polyhedron, and the like.

SUMMARY OF THE INVENTION

The present invention has been developed to solve the aforementionedvarious problems. That is, an object of the present invention is toprovide an image lens distortion correcting method in which (1) nocorrespondence point needs to be given, (2) any special tool or objectis not used, and (3) all points on an image are used, so that aninternal parameter can automatically be obtained with high precision,and an image lens distortion can be corrected with high precision.

According to the present invention, there is provided an image lensdistortion correcting method comprising: an image printing step (A) forprinting an arbitrary image I₁ in a computer; an image pickup step (B)for taking the printed image I₁ with a camera having a lens distortionand obtaining a photographed image I₂ in the computer; a parameterestimating step (C) for obtaining a parameter θ such that the image I₁is artificially distorted with the parameter θ and an obtained distortedimage I₁ ^(ud) agrees with the photographed image I₂; and an imagecorrecting step (D) for using the obtained parameter θ to correct theimage taken by the camera.

According to the aforementioned method of the present invention, firstthe arbitrary image I₁ (Computer graphics, a picture extracted with ascanner, or an image taken by a digital camera) is prepared as acalibration pattern in the computer. Subsequently, the image I₁ isprinted with a printer having no distortion. The printed image I₁ istaken by a camera to be corrected, and the distorted photographed imageI₂ is obtained. When the certain internal parameter θ is given, theimage I₁ can artificially be distorted. When the artificially distortedimage I₁ ^(ud) is compared with the distorted image I₂ and the imagesaccurately equal with each other, the given parameter θ is desirable,and the obtained parameter θ can be used to correct the camera lensdistortion.

According to a preferred embodiment of the present invention, theparameter θ includes a position correction parameter θ^(u) forcorrecting a position, and a distortion correction parameter θ^(d) forcorrecting the distortion.

The position correction parameter θ^(u) is a parameter for conversion toan image I₁ ^(u) with the corrected position from the image I₁. Theconversion parameter is obtained in a least-squares method such that adifference r=I₁(p)−I₁ ^(u)(p+u) between a luminance value I₁(p) of apoint p in the image I₁ and a luminance value I₁ ^(u)(p+u) of a pointp+u in the image I₁ ^(u) corresponding to the point p is minimizedentirely on the image. The obtained parameter is used as the positioncorrection parameter θ^(u).

Moreover, the distortion correction parameter θ^(d) is a parameter forconversion to the image I_(i) ^(u) from the distorted image I₁ ^(ud).The conversion parameter is obtained in the least-squares method suchthat a difference r=I₂(p)−I₁ ^(u)(f(p)) between a luminance value I₂(p)of the point p in the image I₂ and a luminance value I₁ ^(u)(f(p)) of apoint f(p) in the image I₁ ^(u) corresponding to the point p isminimized entirely on the image. The obtained parameter is used as thedistortion correction parameter θ^(d).

Furthermore, in the image correcting step (D), the obtained distortioncorrection parameter θ^(d) is used to correct the image taken by thecamera.

That is, as shown in FIG. 1, a position in the photographed image I₂ inwhich the original image I₁ is projected is not known. Therefore, twostages are performed in order to obtain the distorted image I₁ ^(ud).That is, first the position corrected image I₁ ^(u) is obtained from theimage I₁, subsequently the image I₁ ^(u) is distorted and the distortedimage I₁ ^(ud) is obtained.

First Stage

First, the parameter θ^(u) for conversion to the image I₁ ^(u) from theimage I₁ is obtained. This does not include any internal parameter. Thepoint in I₁ ^(u) corresponding to the point p in the image I₁ deviatesby u and is represented by p+u. Here, u changes with p and θ^(u). Whenthe image I₁ exactly equals with the photographed image I₂, theluminance value I₁ (p) in p has to be equal to the luminance valueI₂(p+u) in p+u. That is, the luminance difference r=I₁(p)−I₂(p+u) ineach point has to be 0. An evaluation function Σr² is obtained bysumming squared differences for all the points, and the parameter θ^(u)is automatically obtained by repeating calculations such that theevaluation function is minimized.

Second Stage

Subsequently, the parameter θ^(d) for conversion to the image I₁ ^(u)from the distorted image I¹ ^(ud) is obtained. This is the internalparameter such as a distortion coefficient. This stage is similar to thefirst stage, but instead of the conversion to the distorted image I₁^(ud) from the image I₁ ^(u), a reverse conversion to the image I₁ ^(u)from the distorted image I₁ ^(ud) is considered.

The point in image I₁ ^(u) corresponding to the point p in thephotographed image I₂ is represented by f(p). Here, f( ) changes with pand θ^(d).

When the photographed image I₂ exactly equals with the image I₁ ^(u),the luminance value I₂(p) in p has to be equal to the luminance value I₁^(u)(f(p)) in f(p). That is, the luminance difference r=I₂(p)−I₁^(u)(f(p)) in each point has to be 0. The evaluation function Σr² isobtained by summing squared differences for all the points, and theparameter θ^(u) is automatically obtained by repeating calculations suchthat the evaluation function is minimized.

Correction

The internal parameter θ^(d) used for obtaining the distorted image I₁^(ud) from the image I₁ can be used to conversely obtain I₁ from I₁^(ud). That is, a processing for correcting the applied distortion canbe performed.

When the image I₂ taken by the camera is subjected to the sameprocessing as described above, the image with the corrected distortioncan be obtained.

In the parameter estimating step (C), it is preferable to alternatelyrepeatedly use the position correction parameter θ^(u) and thedistortion correction parameter θ^(d) and obtain the parameter θ.

Moreover, the parameter θ is a parameter for conversion to the image I₁^(ud) from the image I₁. The conversion parameter is obtained in theleast-squares method such that a difference r=I₁(p)−I₂(d(p+u(p)) betweena luminance value I₁(p) of the point p in the image I₁ and a luminancevalue I₂(d(p+u(p)) of a point d(p+u(p)) in the image I₂ corresponding tothe point p is minimized entirely on the image. The obtained parameteris used as the parameter θ.

According to the method, the parameter θ can be estimated by a singlestep.

Other objects and advantageous characteristics of the present inventionwill be apparent from the following description with reference to theaccompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram schematically showing a method of the presentinvention.

FIG. 2 shows a halftone image displayed on a display in a calibrationpattern used in an embodiment.

FIG. 3 shows the halftone image displayed on the display before andafter correction according to the method of the present invention.

FIG. 4 shows another halftone image displayed on the display before andafter the correction according to the method of the present invention.

FIG. 5 is a schematic diagram of a method of repeating two steps.

FIG. 6 is a schematic diagram of an estimating method with one step.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

A preferred embodiment of the present invention will be describedhereinafter with reference to the drawings. Additionally, commonportions in the respective drawings are denoted with the same referencenumerals, and redundant description is omitted.

Camera calibration and lens distortion correction are importantprocesses in computer vision. In recent years, self calibration has beenstudied. In many researches, to facilitate the calibration, distortionis not considered and respective problems are formulated. Therefore, itis necessary to calibrate an internal parameter of a camera beforehand.

Some pieces of calibration software have been utilized via Internet(e.g., Tsai's method). However, in such ordinary method, it is necessaryto obtain a correspondence between a point on an image and a point on aknown three-dimensional coordinate (a plane, or a cube, a model oranother structure), and subsequently point conversion is estimated.

In this method, the point correspondence must manually be estimated.Therefore, an error is generated in a manual operation and this lacks inreliability. Additionally, the method requires much time and endurance.Therefore, it is almost impossible to measure a change of a distortionparameter during zoom-in/zoom-out of the camera.

One of alternative methods is detection of a marker. This can beperformed by a template matching method to an extend of sub-pixels.However, there is another problem that the marker on the imagecorresponds to a coordinate in a known space. To improve estimationprecision, it is necessary to increase the number of markers, and thisproblem cannot be ignored. Even if this problem can be avoided, thenumber of points for the correspondence is still limited.

In the present invention, there is proposed a new correcting method forcorrecting an image distortion by a lens zoom operation. In the proposedmethod, the correspondence is set, and the estimation is expected to bemore precise than the marker detection. This is because all points ofthe image are used in the method, instead of using some points of themarker. This method is based on an image registration (superposition)method applied to a motion analysis region, and includes the followingthree steps, that is, affine transformation, plane projectiontransformation, and lens distortion restoration.

Three Estimating Steps by Registration

A basic idea is that a correspondence between points is necessary forcalibration, and registration can satisfy this requirement. In theproposed method, an ideal calibration pattern is associated with adistorted image obtained by taking the printed pattern with the camera.Since a size of the pattern printed on a sheet can easily be measured, athree-dimensional coordinate of each pixel in the pattern can easily bedetermined. As other characteristics of the method, any image can beused as the calibration pattern. Moreover, once the parameter isestimated, the parameter can be used as long as the lens zoom operationis unchanged.

The method includes the following three steps. In a first step, thepattern is coarsely transformed into an image represented by an affineparameter and translatory parameter. In a second step, an accurateparameter of plane transformation under through-view projection isestimated by a nonlinear optimizing method. In a third step, adistortion parameter is estimated in order to minimize a residual whichremains after application of the second step.

Affine Transformation by Detected Marker

First, a pattern image I₁ is prepared. The pattern image has three-color(red, green, blue) markers having coordinates m_(r), m_(g), m_(b) inrespective corners. Subsequently, the pattern image I₁ is printed incolor in order to facilitate detection of the marker.

Therefore, an image I₂ of the printed pattern I₁ is taken by a camera tobe calibrated. Coordinates of markers m_(r)′, m_(g)′, m_(b)′ in theimage I₂ are detected by thresholding or template matching.

Here, a parameter for transforming the pattern image I₁ to the patternprojected in the image I₂ is calculated. The transformation isrepresented by six affine parameters θ^(α)=(θ^(α) ₁, . . . θ^(α) ₆)^(T).

While p=(x, y)^(T) is set as one point on the pattern I₁, p+a(p;θ^(α))is set as one point on I₂ corresponding to p. Here, equation (1) of[Equation 1] is established. Subsequently, a linear equation of equation(2) is solved, and parameter θ^(α) is obtained.

Equation 1 $\begin{matrix}{{a( {p;\theta^{a}} )} = {\begin{pmatrix}x & y & 0 & 0 & 1 & 0 \\0 & 0 & x & y & 0 & 1\end{pmatrix}\theta^{a}}} & (1)\end{matrix}$

 m′ _(i) =m _(i)+α(m _(i); θ^(α)),i=r, g, b  (2)

Image Registration; Perspective Fitting

In the first step described above, only three correspondence points areemployed with respect to the affine transformation. Subsequently, planeprojection transformation is utilized, and accurate correspondence amongindividual points is set as image registration.

Model and Formulation

Here, a luminance residual between p_(i) in the image I₁ and p₁+u(p_(i);θ^(u)) in the image I₂ corresponding to p_(i) is minimized. Here, whenthe residual is defined as in equations (3), (4) of [Equation 2], afunction to be minimized is obtained as in equation (5).

Equation 2

r _(i) =I ₁(p _(i))−I₂(p _(i) +u(p; θ^(u)))  (3) $\begin{matrix}\begin{matrix}{{u( {p;\theta^{u}} )} = {\begin{pmatrix}x & y & 0 & 0 & 1 & 0 & x^{2} & {x\quad y} \\0 & 0 & x & y & 0 & 1 & {x\quad y} & y^{2}\end{pmatrix}\theta^{u}}} \\{= {M^{u}\theta^{u}}}\end{matrix} & (4) \\\begin{matrix}{\min\limits_{\theta^{\quad u}}{\sum\limits_{i}{\rho ( r_{i} )}_{1}}} & {{\rho ( r_{i} )} = r_{i}^{2}}\end{matrix} & (5)\end{matrix}$

Here, u( ) denotes displacement of a point on a plane, when a point ofview is changed from a certain point to another point under perspectiveprojection. The displacement is a motion model on the plane under theperspective projection, and is represented by eight parameters θ^(u)=(θ₁^(u), . . . , θ₈ ^(u))^(T) frequently used in motion analysis.

Minimizing Method

The function (5) is minimized by Gauss-Newton method as a famousnonlinear optimization technique and the parameter θ^(u) is estimated.The parameter is updated by equation (6) of [Equation 3].

As an initial value of first six elements of the parameter θ^(u), theaffine parameter θ^(α) obtained in the first stage is used. Last twoelements of the parameter θ^(u) are initialized/set to 0.

A direction δθ^(u) in which the function (5) is decreased is calculatedby equations (7) to (10).

Equation 3

θ^(u)←θ^(u)+δθ^(u)  (6)

δθ^(u)=−(J{tilde over (D)}J ^(T))⁻¹ J{tilde over (D)}r  (7)$\begin{matrix}{J = {{J( \theta^{u} )} = {\frac{\partial r}{\partial\theta^{u}} = \lbrack \frac{\partial r_{i}}{\partial\theta_{j}^{u}} \rbrack}}} & (8) \\{\overset{\_}{D} = {d\quad i\quad a\quad {g\lbrack \frac{\overset{.}{\rho}( r_{i} )}{r_{i}} \rbrack}}} & (9) \\{{\overset{.}{\rho}( r_{i} )} =  \frac{\partial{\rho (r)}}{\partial r} |_{r = r_{i}}} & (10)\end{matrix}$

This is the same as formulation of the least-squares method. That is,the equation is the same as a simultaneous linear equation representedby equation (11) of [Equation 4] for k=1, . . . , 8. A differentiatingchained law is used to convert a partial derivation as represented byequation (12). The calculation of δθ^(u) in the equation (6) is repeateduntil δθ^(u) converges. In each repetition, the parameter estimated inthe previous repetition is used in calculation of u(p). When therepeated calculation stops, the estimated parameter is θ^(u)′,

Equation 4 $\begin{matrix}{{{\sum\limits_{l,i}{\frac{\overset{.}{\rho}( r_{i} )}{r_{i}}\frac{\partial r_{i}}{\partial\theta_{k}^{u}}\frac{\partial r_{i}}{\partial\theta_{l}^{u}}{\delta\theta}_{l}^{u}}} = {- {\sum\limits_{i}{\frac{\overset{.}{\rho}( r_{i} )}{r_{i}}r_{i}\frac{\partial r_{i}}{\partial\theta_{k}^{u}}}}}}{{{{for}\quad k} = 1},\ldots \quad,8.}} & (11) \\{\frac{\partial r}{\partial\theta^{u}} = {{\frac{\partial u}{\partial\theta^{u}}\frac{\partial r}{\partial u}} = {{- M^{u\quad T}}{\nabla{I_{2}( {p + {u(p)}} )}}}}} & (12)\end{matrix}$

Image Registration; Distortion Fitting

In the previous step, except the influence by the lens distortion, theimage registration between the pattern image and the distorted imagepattern ends.

Distortion Model

A relation between a coordinate having no distortion is and coordinatehaving the distortion in a certain image is usually modeled by thefollowing five camera internal parameters, that is, distortionparameters k₁ and k₂, image center coordinate (c_(x), c_(y))^(T) and aratio s_(x) of a pixel width to height. These parameters are set toθ^(d)=(k₁, k₂, c_(x), c_(y), s_(x))^(T).

The distortion is represented by a coordinate system having an origin in(c_(x), c_(y))^(T). On the other hand, the previously used coordinatesystem has an origin in an uppermost left corner. Therefore, anothernotation is introduced as follows.

That is, P_(u)=(x_(u), y_(u))^(T) and p_(d)=(x_(d), u_(d))^(T) are setas points in an image having no distortion and image having thedistortion, respectively, and both images have an origin in the leftupper corner of the image. Moreover, (ξ_(u), η_(u))^(T) and (ξ_(d),η_(d))^(T) are set as points in the image having no distortion and imagehaving the distortion, respectively, in which respective origins are inthe image center (c_(x), c_(y))^(T). These points can be represented byequations (13) to (15) of [Equation 5].

As described above, ξ_(u) is represented to be positive as the functionof ξ_(d) but this does not apply to ξ_(d). In order to obtain ξ_(d) fromξ_(u), equation (16) is repeatedly solved.

Equation 5

(ξ_(u), η_(u))^(T)=(x _(u) , y _(u))^(T)−(c _(x) , c _(y))^(T)  (13)

(ξ_(d), η_(d))^(T)=(ξ_(u), η_(u))^(T)−(κ₁ R ²+κ₂ R ⁴)(ξ_(d),η_(d))^(T)  (14)

(x _(d) , y _(d))^(T)=(s _(x)ξ_(d), η_(d))^(T)+(c _(x) , c_(y))^(T)  (15)

R={square root over (ξ_(d) ²+η_(d) ²)}

$\begin{matrix}{( {\xi_{d\quad k},\eta_{d\quad k}} )^{T} = \frac{( {\xi_{u},\eta_{u}} )^{T}}{1 + {\kappa_{1}R_{k - 1}^{2}} + {\kappa_{2}R_{k - 1}^{4}}}} & (16)\end{matrix}$

 R _(k)={square root over (ξ_(dk) ²+η_(dk) ²)},

R ₀={square root over (ξ_(u) ²+η_(u) ²)}

The aforementioned relation is used to obtain two functions betweenp_(u) and p_(d) for the coordinate system having the origin in theuppermost left corner as represented by equations (17) and (18) of[Equation 6]. Here, f and d are mutually inverse functions. However,since d corresponds to the equation (16), d is not a explicit functionof p_(u).

In any case, the equations (17) and (18) can be used to representtransformation between two images. While I₁ ^(u) is the image of thepattern transformed by applying θ^(u)′ to I₁, I₁ ^(ud) is the image ofthe pattern transformed by applying θ^(d) to I₁ ^(u). That is, therelation is represented by equations (20) to (22). These relations areshown in FIG. 1.

Equation 6

P _(d) =d(p _(u); θ^(d))  (17)

$\begin{matrix}\begin{matrix}{p_{u} = {f( {p_{d};\theta^{d}} )}} \\{= \begin{pmatrix}{{\frac{x_{d} - c_{r}}{s_{r}}( {1 + {\kappa_{1}R^{\prime 2}} + {\kappa_{2}R^{\prime 4}}} )} + c_{x}} \\{{( {y_{d} - c_{y}} )( {1 + {\kappa_{1}R^{\prime 2}} + {\kappa_{2}R^{\prime 4}}} )} + c_{y}}\end{pmatrix}}\end{matrix} & (18) \\{R^{\prime} = \sqrt{( \frac{x_{d} - c_{r}}{s_{r}} )^{2} + \quad ( {y_{d} - c_{y}} )^{2}}} & (19)\end{matrix}$

 I _(i)(p)=I _(i) ^(u)(p+u(p; θ^(u)))  (20)

I ₁ ^(u)(p)=I₁ ^(ud)(d(p))  (21)

I _(i) ^(u)(f(p))=i ₁ ^(ud)(p)  (22)

Minimization by Reverse Registration

The transformation to the image I₁ ^(ud) from I₁ is represented by u andd. However, since d is not the explicit function, the parameter betweenI₁ and I₁ ^(ud) cannot directly be estimated by registration.

Here, reverse registration is considered. As shown in FIG. 1, the imageI₂ is matched with I₁ ^(ud). That is, the residual of luminance betweentwo images represented by equation (23) of [Equation 7] is minimized.

However, this equation can be represented by equation (24) usingequation (22). Then, an estimating method is the same as that of theprevious step. The minimization is performed on the function of equation(25) by the Gauss-Newton method.

Here, Ω={i;p_(i)εI₂, ∃pεI₁, f(p_(i))=p+u(p)} means that the minimizationhas to be applied to the point inside the region corresponding to thepattern I₁ among points in I₂.

Equation 7

r _(i) =I ₁ ^(ud)(p _(i))−I ₂(p _(i))  (23)

r ₁ =I ₂(p _(i))−I ₁ ^(u)(f(p _(i); θ^(d))  (24)

$\begin{matrix}{{\min\limits_{\theta^{d}}{\sum\limits_{\Omega}{\rho ( r_{i} )}}}{\Omega = \{ {{i;{p_{i}\varepsilon \quad I_{2}}},{\exists{p\quad \varepsilon \quad I_{1}}},{{f( p_{i} )} = {p + {u(p)}}}} \}}} & (25)\end{matrix}$

The simultaneous linear equation in which solution needs to be found hasthe same form as that of the equation (11), and is represented byequation (26) of [Equation 8]. The derivation in the equation (26) isrepresented by equation (27).

Equation 8 $\begin{matrix}{{\sum\limits_{l,\quad {i\quad {\varepsilon\Omega}}}{\frac{\overset{.}{\rho}( r_{i} )}{r_{i}}\frac{\partial r_{i}}{\partial\theta_{k}^{d}}\frac{\partial r_{i}}{\partial\theta_{l}^{d}}{\delta\theta}_{l}^{d}}} = {- {\sum\limits_{i\quad {\varepsilon\Omega}}{\frac{\overset{.}{\rho}( r_{i} )}{r_{i}}r_{i}\frac{\partial r_{i}}{\partial\theta_{k}^{d}}}}}} & (26) \\{\frac{\partial r}{\partial\theta^{d}} = {{\frac{\partial f}{\partial\theta^{d}}\frac{\partial r}{\partial f}} = {\frac{\partial f}{\partial\theta^{d}}( {- {\nabla{I_{1}^{u}( {f(p)} )}}} )}}} & (27)\end{matrix}$

Jacobian matrix of the equation (18) is represented by equations (28) to(34) of [Equation 9].

Initial parameters c^(x), c_(y), k₂ and s_(x) for solving the equation(26) are set to ½ of width and height of I₂, 0 and 1, respectively. Onthe other hand, in all of equations (29) to (32), k₁ is arbitrarilyinitially set to avoid 0 by initial setting of k₁=k₂=0. Here, k₁ ∈[−10⁻⁷, 10⁻⁷] is empirically selected.

Equation 9 $\begin{matrix}{\frac{\partial{f( p_{d} )}}{\partial\theta^{d}} = \begin{pmatrix}{{R^{\prime 2}x_{d}} - \frac{c_{r}}{s_{r}}} & {R^{\prime 2}( {y_{d} - c_{y}} )} \\{{R^{\prime 4}x_{d}} - \frac{c_{r}}{s_{r}}} & {R^{\prime 4}( {y_{d} - c_{y}} )} \\\frac{\partial x_{y}}{\partial c_{r}} & \frac{\partial y_{u}}{\partial c_{r}} \\\frac{\partial x_{u}}{\partial c_{y}} & \frac{\partial y_{u}}{\partial c_{y}} \\\frac{\partial x_{u}}{\partial s_{r}} & \frac{\partial y_{u}}{\partial s_{r}}\end{pmatrix}} & (28) \\\begin{matrix}{\frac{\partial x_{u}}{\partial c_{x}} = {1 - {\frac{1}{s_{x}}( {1 + {\kappa_{1}R^{\prime 2}} + {\kappa_{2}R^{\prime 4}}} )}}} \\{{{- 2}( {\kappa_{1} + {2\kappa_{2}R^{\prime 2}}} )\frac{( {x_{d} - c_{x}} )^{2}}{s_{x}^{3}}}}\end{matrix} & (29) \\{\frac{\partial y_{u}}{\partial c_{x}} = {{- 2}( {\kappa_{1} + {2\kappa_{2}R^{\prime 2}}} )\frac{x_{d} - c_{x}}{s_{x}^{2}}( {y_{d} - c_{y}} )}} & (30) \\{\frac{\partial x_{u}}{\partial c_{y}} = {{- 2}( {\kappa_{1} + {2\kappa_{2}R^{\prime 2}}} )\frac{x_{d} - c_{x}}{s_{x}}( {y_{d} - c_{y}} )}} & (31) \\\begin{matrix}{\frac{\partial y_{u}}{\partial c_{y}} = {1 - ( {1 + {\kappa_{1}R^{\prime 2}} + {\kappa_{2}R^{\prime 4}}} )}} \\{{{- 2}( {y_{d} - c_{y}} )^{2}( {\kappa_{1} + {2\kappa_{2}R^{\prime 2}}} )}}\end{matrix} & (32) \\\begin{matrix}{\frac{\partial x_{u}}{\partial s_{x}} = {\frac{- ( {x_{d} - c_{x}} )}{s_{x}^{2}}( {1 + {\kappa_{1}R^{\prime 2}} + {\kappa_{2}R^{\prime 4}}} )}} \\{{{- 2}( {\kappa_{1} + {2\kappa_{2}R^{\prime 2}}} )\frac{( {x_{d} - c_{x}} )^{3}}{s_{x}^{4}}}}\end{matrix} & (33) \\{\frac{\partial y_{u}}{\partial s_{x}} = {{- 2}( {y_{d} - c_{y}} )( {\kappa_{1} + {2\kappa_{2}R^{\prime 2}}} )\frac{( {x_{d} - c_{x}} )^{2}}{s_{x}^{3}}}} & (34)\end{matrix}$

Repetition of Two Steps

FIG. 5 is a schematic diagram showing a method of repeating two steps.Even when there is a relatively large residual in an initial state, astrategy for shifting to a fine resolution from a coarse resolution(coarse-to-fine) is employed in order to shorten a calculation time andperform precise estimation. That is, first a filtered and blurred imageis processed, a filter is gradually changed to a filter having a smallerblur, and a cleared image is processed. Therefore, while the resolutionsof I₁, I₁ ^(u), I₂ are changed to the fine resolution from the coarseresolution, second and third steps are alternately repeated (see FIG.5).

Therefore, the estimation of the second step has to be corrected toinclude an estimated value θ^(d) of the third step.

Here, I₂ ^(f) is assumed to be an image obtained by correctingdistortion of the photographed image I₂ through transformation f usingthe parameter θ^(d). That is, the image can be represented by equations(35), (36).

This is used to correct the equation (3) of the residual as representedby equation (37).

In this case, Jacobian equation (12) turns to equation (38).

The corrected second step, and the third step are alternately repeateduntil the estimated value converges.

The image having no distortion can be obtained from the photographedimage by the transformation f using the parameter θ^(d) out of theobtained parameters.

Equation 10

I ₂ ^(ƒ)(ƒ(p ^(; θd)); θ^(d))=I ₂(p)  (35)

I ₂ ^(f)(p ^(; θd))=I ₂(d(p ^(; θd)))  (36)

r _(i) =I ₁(p _(i))−I ₂ ^(ƒ)(p _(i) +u(p _(i) ; θu); θ^(d))  (37)

$\begin{matrix}{\frac{\partial r}{\partial\theta^{u}} = {{- M^{u\quad T}}{\nabla{I_{2}^{f}( {{p + {u( {p;\theta^{u}} )}};\theta^{d}} )}}}} & (38)\end{matrix}$

Estimation with One Step

FIG. 6 is a schematic diagram showing an estimating method with onestep. The estimation by the repetition of two steps has been describedabove, but these two steps can be unified into one estimating step.

Since d is not the explicit function, that is, d is an implicit functionunable to be represented by any equation, a function obtained bysynthesizing u and d also becomes implicit function. This synthesizedfunction u·d transforms the image I₁ to I₁ ^(ud). Therefore, theparameter θ^(u), θ^(d) may be obtained such that I^(ud) accuratelyequals with I₂. In this case, Jacobian of the synthesized function needsto be explicit, and the explicit function is obtained from Jacobian of uand d.

In actual, since d is an implicit function, Jacobian cannot be obtainedby directly differentiating the equation. However, when the implicitfunction theorem known in a mathematical field is used, Jacobian of dcan be obtained as the explicit function.

Here, it is assumed that p_(u) and θ^(d) are unified into one, andrepresented by q=(p_(u), θ^(d)) Then, function F of equation (39) as thefunction of q is considered.

In this case, d is an implicit function defined by p_(d)=d(q).Therefore, Jacobian of d can be calculated by the implicit functiontheorem as represented by equations (40) and (41).

Here, E₂ is a 2×2 unit matrix, and represented by equation (41a).Moreover, ∂f/∂θ^(d) is the same as the equation (28). Moreover,∂f/∂p_(d) can be calculated from the equation (18), and is actuallyrepresented by equations (42) and (43) using the equations (29) to (32).

Equation 11

F(q,p _(d))=p _(u)−ƒ(p _(d), θ^(d))=0  (39)

$\begin{matrix}{\frac{\partial d}{\partial q} = {{- \frac{\partial F^{- 1}}{\partial p_{d}}}\frac{\partial F}{\partial q}}} & (40) \\{\quad {= {{- \frac{\partial f^{- 1}}{\partial p_{d}}}( {{- E_{2}}\frac{\partial f}{\partial\theta^{d}}} )}}} & (41) \\{E_{2} = \begin{pmatrix}1 & 0 \\0 & 1\end{pmatrix}} & ( \text{41a} ) \\{\frac{\partial f}{\partial p_{d}} = \begin{pmatrix}\frac{\partial x_{u}}{\partial x_{d}} & \frac{\partial y_{u}}{\partial x_{d}} \\\frac{\partial x_{u}}{\partial y_{d}} & \frac{\partial y_{u}}{\partial y_{d}}\end{pmatrix}} & (42) \\{\quad {= \begin{pmatrix}{1 - \frac{\partial x_{u}}{\partial c_{x}}} & {- \frac{\partial y_{u}}{\partial c_{x}}} \\{- \frac{\partial x_{u}}{\partial c_{y}}} & {1 - \frac{\partial y_{u}}{\partial c_{y}}}\end{pmatrix}}} & (43)\end{matrix}$

As described above, Jacobian of the function obtained by synthesizingtransformations u, d can be calculated. A squared sum of luminancedifferences from the point in the photographed image I₂ corresponding tothe point p_(i) in the pattern I₁, and is represented by equations (44),(45).

Parameter θ={θ^(u), θ^(d)} for minimizing the squared sum is estimated(see FIG. 6). The optimizing method for the estimation is performedsimilarly as the estimation of the second and third steps. That is, thesimultaneous linear equation to be solved is represented by equation(46).

Here, partial differential ∂r_(i)/∂θ_(k) is derived as an element ofJacobian of equations (47) to (52).

Here, p_(d)=d(p_(u)), p_(u)=p+u (p) Moreover, E₅ is a 5×5 unit matrix.

The image having no distortion can be obtained from the photographedimage by the transformation f using the parameter θ^(d) out of theobtained parameters.

Equation 12

r _(i) =I ₁(p _(i))−I ₂(d(p _(i) +u(p _(i) ^(, θu))), θ^(d))  (44)

$\begin{matrix}{{\min\limits_{\theta^{u}}{\sum\limits_{i}{\rho ( r_{i} )}}},{{\rho ( r_{i} )} = r_{i}^{2}}} & (45) \\{{\sum\limits_{l}{\sum\limits_{i}{\frac{\overset{.}{\rho}( r_{i} )}{r_{i}}\frac{\partial r_{i}}{\partial\theta_{k}}\frac{\partial r_{i}}{\partial\theta_{l}}{\delta\theta}_{l}}}} = {- {\sum\limits_{i}{\frac{\overset{.}{\rho}( r_{i} )}{r_{i}}r_{i}\frac{\partial r_{i}}{\partial\theta_{k}}}}}} & (46) \\{\frac{\partial r}{\partial\theta} = {\frac{\partial r}{\partial q}\frac{\partial q}{\partial\theta}}} & (47) \\{= {\frac{\partial r}{\partial d}\frac{\partial d}{\partial q}( {\frac{\partial q}{\partial\theta^{u}}\frac{\partial q}{\partial\theta^{d}}} )}} & (48) \\{= {\frac{\partial r}{\partial d}\frac{\partial d}{\partial q}\begin{pmatrix}\frac{\partial p_{u}}{\partial\theta^{u}} & \frac{\partial p_{u}}{\partial\theta^{d}} \\\frac{\partial\theta^{d}}{\partial\theta^{u}} & \frac{\partial\theta^{d}}{\partial\theta^{d}}\end{pmatrix}}} & (49) \\{= {\frac{\partial r}{\partial d}\frac{\partial d}{\partial q}\begin{pmatrix}\frac{\partial p_{u}}{\partial\theta^{u}} & 0 \\0 & E_{5}\end{pmatrix}}} & (50) \\{= {{- {\nabla\quad {I_{2}( {d( p_{u} )} )}}}\frac{\partial d}{\partial q}\begin{pmatrix}M^{u} & 0 \\0 & E_{5}\end{pmatrix}}} & (51) \\{= {{\nabla\quad {I_{2}( {d( p_{u} )} )}}\frac{\partial f^{- 1}}{\partial p_{d}}\begin{pmatrix}{- M^{u}} & \frac{\partial f}{\partial\theta^{d}}\end{pmatrix}}} & (52)\end{matrix}$

Several Strategies

(1) Interpolation of Pixel Value

When the luminance of the point having no coordinate on an integer gridis necessary, bilinear interpolation is applied between adjacent pixelvalues.

(2) Histogram Matching

For the image taken by the camera, since pattern luminance oftenchanges, image histogram is converted so that the histogram of the imagebecomes the same as that of the pattern.

EXAMPLES

While the camera zoom parameter was changed, an experiment was performedfor the actual image taken by the camera in the proposed method. Apicture shown in FIG. 2 was used as the calibration pattern printed by alaser printer (EPSON LP-9200PS2), and an image of the picture was takenby a CCD camera (Sony EVI-D30) fixed on a support with image takingsoftware (SGI O2). While the camera zoom was changed from a widestangle, two images of grid and calibration patterns were photographed.The wider the angle is, the more the distortion influence increases.Then, a grid line is curved.

An estimation result of the internal parameter according to the presentinvention is shown in Table 1. Moreover, this internal parameter is usedto show images before and after the correction in FIG. 3 and FIG. 4. Itis seen that the curved line in the grid pattern is transformed into astraight line by the correction.

TABLE 1 FIG. 3 (a) FIG. 3 (e) k₁ 2.804e-07 −6.7631e-08 k₂ 2.992e-135.219e-13 c_(x) 327.8 326.6 c_(y) 214.3 184.2 s_(x) 0.9954 0.9997

A gradation change of an image illuminance cannot be removed by a simplehistogram conversion which is to be replaced with the estimation of theilluminance change by some methods such as linear luminance restriction.Therefore, the grid line is still slightly curved particularly in aperipheral portion of the image, and as shown on the right side of FIG.3, there is a slight error in the experiment result.

In an experiment simulation in which a noise was added to each pixel ofthe transformed pattern and the pattern was used as the distorted image,the proposed method satisfactorily functioned even with a uniformamplitude of added noise of ±50 or more.

As described above, in the present invention, there is proposed aninventive method for obtaining the internal parameter of the camera tocorrect the image distortion. The proposed method is based on the imageregistration, and includes two nonlinear optimizing steps of perspectivefitting and distortion fitting by geometric transformation. Theexperiment result has proved an efficiency of the present proposedmethod which can reduce an operator's labor. The nonlinear optimizationrequires time, but the technique is sufficiently effective as a batchprocess.

As described above, there is proposed an inventive image lens distortioncorrecting method for obtaining the camera internal parameter to correctthe distortion of the image taken by the camera with the zoom lensattached thereto. This method is based on image registration. First, acertain calibration pattern is transformed into the distorted imagepattern using affine transformation. Subsequently, the nonlinearoptimization of Gauss-Newton method for minimizing the luminancedifference of pixel values of two images is applied, and the two imagesare registered. Finally, the distortion parameter is estimated tominimize the luminance difference remaining after the second step. Theexperiment result shows effectiveness of the proposed method.

Therefore, in the image lens distortion correcting method of the presentinvention, (1) no correspondence point needs to be given, (2) anyspecial tool or object is not used, and (3) all points on the image areused, so that the internal parameter can automatically be obtained withhigh precision, the image lens distortion can be corrected with highprecision, and other superior effects are produced.

Additionally, the present invention is not limited to the aforementionedembodiment, and can of course variously be changed within the scope ofthe present invention.

What is claimed is:
 1. An image lens distortion correcting methodcomprising: an image printing step (A) for printing a preparedcalibration pattern image I₁ in a computer; an image pickup step (B) forphotographing the printed image I₁ with a camera having a lensdistortion and obtaining a photographed image I₂ in the computer; aparameter estimating step (C) for obtaining a parameter θ such that theimage I₁ is artificially distorted with the parameter θ and so as toobtain a distorted image I₁ ^(ud) equal with the photographed image I₂;and an image correcting step (D) for using the parameter θ to correctthe photographed image I₂.
 2. The image lens distortion correctingmethod according to claim 1 wherein said parameter θ comprises aposition correction parameter θ^(u) for correcting a position, and adistortion correction parameter θ^(d) for correcting the distortion. 3.The image lens distortion correcting method according to claim 2,wherein said position correction parameter θ^(u) is a first parameterfor conversion to an image I₁ ^(u) with the corrected position from theimage I₁, the first conversion parameter is calculate in a least-squaresmethod such that a difference r=I₁(p)−I₁ ^(u)(p+u) between a luminancevalue I_(i)(p) of a point p in the image I₁ and a luminance value I₁^(u)(p+u) of a point p+u in the image I₁ ^(u) corresponding to the pointp is minimized entirely on the image I₁, and the first parameter thuscalculated is used as the position correction parameter θ^(u).
 4. Theimage lens distortion correcting method according to claim 3, whereinsaid distortion correction parameter θ^(d) is a second parameter forconversion to the image I₁ ^(u) from the distorted image I₁ ^(ud), thesecond conversion parameter is calculated in the least-squares methodsuch that a difference r=I₂(p)−I₁ ^(u)(f(p)) between a luminance valueI₂(p) of the point p in the image I₂ and a luminance value I₁ ^(u)(f(p))of a point f(p) in the image I₁ ^(u) corresponding to the point p isminimized entirely on the image I₂, and the second parameter thuscalculated is used as the distortion correction parameter θ^(d).
 5. Theimage lens distortion correcting method according to claim 4, whereinsaid parameter estimating step (C) comprises a step of alternatelyrepeatedly using said position correction parameter θ_(u) and saiddistortion correction parameter θ^(d) to calculate the parameter θ. 6.The image lens distortion correcting method according to claim 1,wherein said parameter θ is a parameter for conversion to the image I₁^(ud) from the image I₁, the conversion parameter is calculated in aleast-squares method by using an implicit function theorem such that adifference r=I₁(p)−I₂(d(p+u(p)) between a luminance value I₁(p) of apoint p in the image I₁ and a luminance value I₂(d(p+u(p)) of a pointd(p+u(p)) in the image I₂ corresponding to the point p is minimizedentirely on the image I₁, and the conversion parameter thus calculatedis used as the parameter θ.
 7. The image lens distortion correctingmethod according to claim 4, wherein said image correcting step (D)comprises a step of using the calculated distortion correction parameterθ^(d) to correct the photographed image I₂.
 8. An image lens distortioncorrecting method comprising the steps of: (a) printing a preparedcalibration pattern image I₁ in a computer; (b) photographing theprinted image I₁ with a camera having a lens distortion to obtain adistorted photographed image I₂ in the computer; (c) calculating aparameter θ such that the image I₁ is artificially distorted with theparameter θ so as to obtain a distorted image I₁ ^(ud) equal with thephotographed image I₂; and (d) using the calculated parameter θ tocorrect the photographed image I₂.
 9. An image lens distortioncorrecting method according to claim 8, wherein said parameter θcomprises a position correction parameter θ^(u) for correcting aposition, and a distortion correction parameter θ_(d) for correctingdistortion.
 10. An image lens distortion correcting method according toclaim 9, wherein said position correction parameter θ^(u) is a firstparameter for conversion to an image I₁ ^(u) with the corrected positionfrom the image I₁, the first conversion parameter is calculated in aleast-squares method such that a difference r=I₁(p)−I₁ ^(u)(p+u) betweena luminance value I₁(p) of a point p in the image I₁ and a luminancevalue I₁ ^(u)(p+u) of a point p+u in the image I₁ ^(u), corresponding tothe point p, is minimized entirely on the image I₁, and the firstconversion parameter thus calculated is used as the position correctionparameter θ^(u).
 11. An image lens distortion correcting methodaccording to claim 10, wherein said distortion correction parameterθ^(d) is a second parameter for conversion to the image I₁ ^(u) from thedistorted image I₁ ^(ud), the second conversion parameter is calculatedin the least-squares method such that a difference r=I₂(p)−I₁ ^(u)(f(p))between a luminance value I₂(p) of the point p in the image I₂ and aluminance value I₁ ^(u)(f(p)) of a point f(p) in the image I₁ ^(u),corresponding to the point p, is minimized entirely on the image I₂, andthe second conversion parameter thus calculated is used as thedistortion correction parameter θ^(d).
 12. An image lens distortioncorrecting method according to claim 11, wherein calculation of theparameter θ comprises alternately repeatedly using said positioncorrection parameter θ^(u) and said distortion correction parameterθ^(d) to calculate the parameter θ.
 13. An image lens distortioncorrecting method according to claim 8, wherein said parameter θ is aparameter for conversion to the image I₁ ^(ud) from the image I₁ that iscalculated in a least-squares method by using an implicit functiontheorem such that a difference r=I₁(p)−₂(d(p+u(p)) between a luminancevalue I₁(p) of a point p in the image I₁ and a luminance valueI₂(d(p+u(p)) of a point d(p+u(p)) in the image I₂, corresponding to thepoint p, is minimized entirely on the image I₁, and the parameter forconversion to the image I₁ ^(ud) from the image I₁ thus calculated isused as the parameter θ.
 14. The image lens distortion correcting methodaccording to claim 11, wherein the calculated distortion correctionparameter θ^(d) is used to correct the photographed image I₂.
 15. Theimage lens distortion correcting method according to claim 9, whereinsaid distortion correction parameter θ^(d) is a five parameterdistortion comprising k₁, k₂, c_(x), c_(y), and s_(x), where k₁ and k₂are distortion parameters, c_(x) and c_(y) represent an image centercoordinate, and s_(x) is a ratio of pixel width to pixel height.
 16. Theimage lens distortion correcting method according to claim 15, whereinsaid position correction parameter θ^(u) is an eight parameterdisplacement comprising θ₁ ^(u), θ₂ ^(u), θ₃ ^(u), θ₄ ^(u), θ₅ ^(u), θ₆^(u), θ₇ ^(u), and θ₈ ^(u).